close all
clear
clc
%%%%%%%%%%%%%%%%%%%%%
%     Parameter     %
%%%%%%%%%%%%%%%%%%%%%
a=0.2;                      %lowpass filter
slideWindow=5;             %delta function forwarding and tailing window
cutoff=20;                  %peak selection cutoff critias
excludePercentage=0.1;      %fitting excluding percentage
fitWindow=5;
cutoffForce=11.01;
stepCutoffMagnet=8.5;
internalTTest=1;            %Single T test
internalTTestAlpha=0.01;
position_overstretching=12.95;
aa_length=0.5;    % the contour length of an individual amino acid
filelist=importdata('G:\RecA filament\ss-recafilament-interaction\filelist.txt');

for currentfile = filelist'

name=char(currentfile);


%% Data processing
steps_data=[];
filename=sprintf('%s%s',name,'.txt')
c=65/(exp(-(13.5-position_overstretching)/0.36)+0.48*exp(-(13.5-position_overstretching)/1.12));
%c=280;
datastruct=importdata(filename);
dataraw=datastruct.data;
egxT=dataraw(:,1)';
egMagneticT=dataraw(:,2);
egRawDynT=dataraw(:,3);
%plot(egxT); % Check equal spacing

%Mingxi data: A. linear drift remover
%             B. force altering regions
%             C. points for fitting selection
linearFitP=polyfit(egxT(egMagneticT==min(egMagneticT))',egRawDynT(egMagneticT==min(egMagneticT)),1);
egfittedDriftT=polyval(linearFitP, egxT)';
egRawDynTShift=egRawDynT-egfittedDriftT;
%figure;
%plot(egxT,egRawDynT);
%hold on; plot(egxT,egfittedDriftT,'r','LineWidth',3);
%figure;
%plot(egxT,egRawDynT-egfittedDriftT,'b');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%cutoffForce=min(egMagneticT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tsig=egMagneticT>cutoffForce;
dsig=diff([0;tsig;0]);
startIndex=find(dsig > 0);
endIndex=find(dsig < 0)-1;


for forceAlteringRegionIndex=1:length(startIndex)
egx=egxT(startIndex(forceAlteringRegionIndex):endIndex(forceAlteringRegionIndex));
egRawDyn=egRawDynTShift(startIndex(forceAlteringRegionIndex):endIndex(forceAlteringRegionIndex));
%scaled magnetic Location
egMagnetic=egMagneticT(startIndex(forceAlteringRegionIndex):endIndex(forceAlteringRegionIndex));
egTruDyn=(egMagnetic-min(egMagnetic))*(max(egRawDyn)-min(egRawDyn))/(max(egMagnetic)-min(egMagnetic));
%figure;
%plot(egxT,egRawDynTShift);
%hold on;plot(egx,egRawDyn,'r');
        

%% Assign
x=egx;
RawDyn=egRawDyn;
%TruDyn=egTruDyn;
%% Filtering By lowpass filter (Freqency**)
% Without filter seems better for loading rate data (or change parameter)
 FilDynA=filter(a,[1 a-1],RawDyn);
 FilDynB=smooth(RawDyn,slideWindow,'sgolay'); % SG smooth
% 
 hFig = figure();
 set(hFig, 'Visible', 'off')
 subplot(3,1,1);
 plot(x,RawDyn,'b',x,FilDynA,'r');
 hold on; plot(x,FilDynB,'g');
 

%% auxiliary function (delta function)

disp([num2str(slideWindow), ' points delta function ...']);
customDyn=FilDynB;
%customDyn=RawDyn;
deltaF=zeros(1,numel(customDyn));

%starting
for i=1:slideWindow
    diff_value=i-1;
    forwarding=mean(customDyn(i:i+diff_value));
    tailing=mean(customDyn(i-diff_value:i));
    deltaF(i)=tailing-forwarding;
end

for i=slideWindow+1:numel(customDyn)-slideWindow
    %disp(i);
    forwarding=mean(customDyn(i:i+slideWindow));
    tailing=mean(customDyn(i-slideWindow:i));
    deltaF(i)=tailing-forwarding;
end

%ending
for i=numel(customDyn)-slideWindow+1:numel(customDyn)
    diff_value=numel(customDyn)-i;
    forwarding=mean(customDyn(i:i+diff_value));
    tailing=mean(customDyn(i-diff_value:i));
    deltaF(i)=tailing-forwarding;
end


%% peak location determination

disp([num2str(cutoff), ' cutoff peak locating ...']);

deltaFCut=deltaF.*(deltaF>=cutoff|deltaF<=-cutoff);
deltaFCutU=deltaF.*(deltaF>=cutoff);  %positive
deltaFCutD=deltaF.*(deltaF<=-cutoff); %negative


[Upeaks,Ulocs]=findpeaks(deltaFCutU,'minpeakdistance',slideWindow);
[Dpeaks,Dlocs]=findpeaks(-deltaFCutD,'minpeakdistance',slideWindow);

subplot(3,1,3);
plot(x,deltaF,'k-');
hold on;plot(x,deltaFCut,'g');
plot(Ulocs*(max(x)-min(x))/(numel(x)-1)+min(x),Upeaks,'ro');
plot(Dlocs*(max(x)-min(x))/(numel(x)-1)+min(x),-Dpeaks,'mo');

%%%%%%%%%%%%%%%%%%%%%%
%        TRAIL       %
%%%%%%%%%%%%%%%%%%%%%%
%% T Test filter remove false steps
fitDyn=RawDyn;
locs=sort([1,Ulocs, Dlocs, numel(fitDyn)]);

% Multiple rounds of test until fix
if internalTTest == 1
    compareDyn=fitDyn;
    compareDyn_size=size(compareDyn);
    templocs=zeros(1,numel(locs));
    previouslocs=0;
    count=1;
    
    while numel(previouslocs) ~= numel(locs)
        disp([num2str(internalTTestAlpha*100), '% Internal Student T Test: ', num2str(count)]);
        previouslocs = locs;
        for i=2:numel(locs)-1
            %interval selection
            startI=locs(i)-fitWindow;
            endI=locs(i)+fitWindow;
            if startI < 1
                startI=1;
            end
            if endI > compareDyn_size
                endI=compareDyn_size;
            end
                
            jumpI=locs(i);
            startCutI=startI+floor((jumpI-startI)*excludePercentage);
            endCutI=endI-floor((endI-jumpI)*excludePercentage);
            groupA=compareDyn(startCutI:jumpI-1);
            groupB=compareDyn(jumpI+1:endCutI);

            [H,PI]=ttest2(groupA,groupB,internalTTestAlpha,[],'unequal');
            %disp(PI);

            if H==1
                templocs(i)=locs(i);
            else
                locs(i)=startI;
            end
        end

        locs=[1,templocs(templocs~=0), numel(fitDyn)];
        templocs=zeros(1,numel(locs));
        count=count+1;
    end
end        

%% least square fit to heaviside function (per interval)
%Change to fit locally to adjacentpoitns
fittedDyn=zeros(numel(fitDyn),1);

for i=2:numel(locs)-1
    %tempfix for some locs<0 error
    if locs(i)>fitWindow && locs(i)<(length(x)-fitWindow)
        %interval selection
        startI=locs(i)-fitWindow;
        endI=locs(i)+fitWindow;
        jumpI=locs(i);
        startCutI=startI+floor((jumpI-startI)*excludePercentage);
        endCutI=endI-floor((endI-jumpI)*excludePercentage);
        currentX=x(startCutI:endCutI);
        currentY=fitDyn(startCutI:endCutI);
        
        %figure;plot(x(startI:endI),fitDyn(startI:endI),'go',currentX,currentY,'ro')
        %guess starting points
        ga=mean(fitDyn(jumpI:endCutI))-mean(fitDyn(startCutI:jumpI));
        gb=x(jumpI);
        if ga>=0
            gc=mean(fitDyn(startCutI:jumpI));
        else
            gc=mean(fitDyn(jumpI:endCutI));
        end
        
        %fit
        fequation='a.*heaviside(x-b)+c';
        ft = fittype(fequation);
        fo = fitoptions('Method','NonlinearLeastSquares','StartPoint',[ga gb gc]);
        currentFit = fit(currentX',currentY,ft,fo);
        
        %currentF=currentFit(currentX);
        %hold on;plot(currentX,currentF,'k','LineWidth',3);
        %currentFit
        
        %update result
        locs(i)=round((currentFit.b-min(x))/(max(x)-min(x))*(numel(x)-1));
        
        leftValue=currentFit.c;
        rightValue=currentFit.a+currentFit.c;
        stepValue=0.5*currentFit.a+currentFit.c;
        if locs(i)>fitWindow
            fittedDyn((locs(i)-fitWindow):locs(i)-1)=leftValue;
        else
            fittedDyn(1:locs(i)-1)=leftValue;
        end
        fittedDyn((locs(i)+1):locs(i)+fitWindow)=rightValue;
        fittedDyn(locs(i))=stepValue;
        %get the parameters for the magnet and calculated the force
        stepsize=currentFit.a;
        locmagnet=egMagnetic(locs(i));
        force_then=c*(exp(-(13.5-locmagnet)/0.36)+0.48*exp(-(13.5-locmagnet)/1.12));
        full_length=stepsize/zz_wlc(force_then*aa_length/4);
        %record data for forces >9.5 to reduce false positives
        if locmagnet>stepCutoffMagnet
            steps_data=[steps_data;[forceAlteringRegionIndex,locmagnet,force_then,stepsize,full_length]];
        end
    end
    
end



% %update last step & first value
% rightValue=currentFit.c+currentFit.a;
% fittedDyn(locs(i)+1:numel(fittedDyn))=rightValue;
% fittedDyn(1)=fittedDyn(2);

% %update locs value
% for i=2:numel(locs)-1
%     fittedDyn(locs(i))=.5*(fittedDyn(locs(i)-1)+fittedDyn(locs(i)+1));
% end

subplot(3,1,2);
plot(x,fitDyn,'ro');
hold on; plot(x,fittedDyn,'g','LineWidth',3);
end
figure_filename = sprintf('%s%s',name,'.png');

print(gcf, '-dpng', figure_filename, '-r 640');
%open write file
filename = sprintf('%s%s',name,'-f-step.dat');
fid = fopen(filename,'w');
if(fid == 0)
    disp('file cannot be opened');
end
fprintf(fid,'%s\n',name);

fprintf(fid,'rampnum\tmagnet(mm)\tforce (pN)\tstepsize(nm)\tconter length(nm)\n');
fprintf(fid,'%d\t%.3f\t%.3f\t%.3f\t%.3f\n',steps_data');

end